rm(list=ls())
library(tidyverse)
library(AER)

##### LOAD MY_experiment DATASET #####
# Set working directory to where datasets are stored
MY_experiment <- read.csv("MY_experiment.csv")

##### RECODE ETHNIC DISCRIMINATION + PERSPECTIVE-TAKING VARIABLES FOR ANALYSIS ##### 
MY_experiment <- MY_experiment %>%
  rowwise() %>%
  mutate(EthnicDis = mean(c(EthnicPM,
                            EthnicParty,
                            EthnicRights,
                            ICERDrat)),
         PT = mean(c(PT1,PT2,PT3,PT4,PT5,PT6,PT7)))

colnames(MY_experiment)

##### DESCRIPTIVE STATISTICS (Table OA3c.29) #####
MY_experiment <- cbind.data.frame(MY_experiment,model.matrix(~Treat2 - 1, data = MY_experiment))
sum_stats <- MY_experiment %>%
  dplyr::select(Treat1,Treat2Bilingual,
                Treat2English,Treat2Malay,
                EthnicDis,PT,
                Female,Education,
                Age,Income,Check,
                EngFluent,
                MalayFluent
                #Duration..in.seconds.
  ) %>%
  psych::describe() %>%
  data.frame() # summary statistics
sum_stats <- sum_stats[,c(2:5, 8:9,11)]

##### BALANCE CHECKS ##### 
###### Balance Checks, Bilingual versus Monolingual (Table OA3d.30) ###### 
bal.table.Treat1 <- MY_experiment %>% 
  group_by(Treat1) %>%
  dplyr::summarize_at(c("Female", "Age", "Education", "Income", "Check",
                        "EngFluent", "MalayFluent"),
                      list(mean = mean, sd = sd),
                      na.rm = T) 
bal.table.Treat1_out <- as.data.frame(t(bal.table.Treat1))
rownames(bal.table.Treat1_out) <- colnames(bal.table.Treat1)
colnames(bal.table.Treat1_out) <- rownames(bal.table.Treat1)
bal.table.Treat1_out <- bal.table.Treat1_out[-1,] 
bal.table.Treat1_out <- bal.table.Treat1_out[order(row.names(bal.table.Treat1_out)), ]
bal.table.Treat1_out$Covariates <- c("Age", "",
                                     "Check", "",
                                     "Education", "",
                                     "Fluent English", "",
                                     "Female", "",
                                     "Income", "",
                                     "Fluent Bahasa Melayu", "")
bal.table.Treat1_out <- dplyr::select(bal.table.Treat1_out,
                                      c(Covariates, "1", "2"))
rownames(bal.table.Treat1_out) <- NULL

b1.Treat1 <- summary(lm(Age ~ Treat1, data = MY_experiment))
b2.Treat1 <- summary(lm(Check ~ Treat1, data = MY_experiment))
b3.Treat1 <- summary(lm(Education ~ Treat1, data = MY_experiment))
b4.Treat1 <- summary(lm(EngFluent ~ Treat1, data = MY_experiment))
b5.Treat1 <- summary(lm(Female ~ Treat1, data = MY_experiment))
b6.Treat1 <- summary(lm(Income ~ Treat1, data = MY_experiment))
b7.Treat1 <- summary(lm(MalayFluent ~ Treat1, data = MY_experiment))

b_fstat.Treat1 <- sapply(list(b1.Treat1, b2.Treat1, b3.Treat1, b4.Treat1,
                              b5.Treat1, b6.Treat1, b7.Treat1), function(x) x$fstatistic[1])
b_fstat_tojoin.Treat1 <- c(sapply(b_fstat.Treat1, c, NA))

b_fstat_pval.Treat1 <- sapply(list(b1.Treat1, b2.Treat1, b3.Treat1, b4.Treat1,
                                   b5.Treat1, b6.Treat1, b7.Treat1), function(x) pf(x$fstatistic[1], x$fstatistic[2], 
                                                                                    x$fstatistic[3], lower.tail = F))
b_fstat_pval_tojoin.Treat1 <- c(sapply(b_fstat_pval.Treat1, c, NA))

bal.table.Treat1_out$b_fstat <- b_fstat_tojoin.Treat1
bal.table.Treat1_out$b_fstat_pval <- b_fstat_pval_tojoin.Treat1

for (i in 2:ncol(bal.table.Treat1_out)){
  bal.table.Treat1_out[,i] <- as.character(round(as.numeric(bal.table.Treat1_out[,i]),3))
  if (i <= 6){
    for (j in 1:nrow(bal.table.Treat1_out)){
      if (j %% 2 == 0){
        bal.table.Treat1_out[j,i] <- paste0("(",bal.table.Treat1_out[j,i],")")
      }
    }
  }
}
colnames(bal.table.Treat1_out) <- c("", "Monolingual", "Bilingual",
                                    "F Statistic", "Prob > F")

###### Balance Checks, Three Experimental Conditions (Table OA3d.31) ###### 
bal.table.Treat2 <- MY_experiment %>% 
  group_by(Treat2) %>%
  dplyr::summarize_at(c("Female", "Age", "Education", "Income", "Check", 
                        "EngFluent", "MalayFluent"),
                      list(mean = mean, sd = sd),
                      na.rm = T) 
bal.table.Treat2_out <- as.data.frame(t(bal.table.Treat2))
rownames(bal.table.Treat2_out) <- colnames(bal.table.Treat2)
colnames(bal.table.Treat2_out) <- rownames(bal.table.Treat2)
bal.table.Treat2_out <- bal.table.Treat2_out[-1,] 
bal.table.Treat2_out <- bal.table.Treat2_out[order(row.names(bal.table.Treat2_out)), ]
bal.table.Treat2_out$Covariates <- c("Age", "",
                                     "Check", "",
                                     "Education", "",
                                     "Fluent English", "",
                                     "Female", "",
                                     "Income", "",
                                     "Fluent Bahasa Melayu", "")
bal.table.Treat2_out <- dplyr::select(bal.table.Treat2_out,
                                      c(Covariates, "1", "2", "3"))
rownames(bal.table.Treat2_out) <- NULL

b1.Treat2 <- summary(lm(Age ~ Treat2, data = MY_experiment))
b2.Treat2 <- summary(lm(Check ~ Treat2, data = MY_experiment))
b3.Treat2 <- summary(lm(Education ~ Treat2, data = MY_experiment))
b4.Treat2 <- summary(lm(EngFluent ~ Treat2, data = MY_experiment))
b5.Treat2 <- summary(lm(Female ~ Treat2, data = MY_experiment))
b6.Treat2 <- summary(lm(Income ~ Treat2, data = MY_experiment))
b7.Treat2 <- summary(lm(MalayFluent ~ Treat2, data = MY_experiment))

b_fstat.Treat2 <- sapply(list(b1.Treat2, b2.Treat2, b3.Treat2, b4.Treat2,
                              b5.Treat2, b6.Treat2, b7.Treat2),
                         function(x) x$fstatistic[1])
b_fstat_tojoin.Treat2 <- c(sapply(b_fstat.Treat2, c, NA))

b_fstat_pval.Treat2 <- sapply(list(b1.Treat2, b2.Treat2, b3.Treat2, b4.Treat2,
                                   b5.Treat2, b6.Treat2, b7.Treat2),
                              function(x) pf(x$fstatistic[1],
                                             x$fstatistic[2], 
                                             x$fstatistic[3], lower.tail = F))
b_fstat_pval_tojoin.Treat2 <- c(sapply(b_fstat_pval.Treat2, c, NA))

bal.table.Treat2_out$b_fstat <- b_fstat_tojoin.Treat2
bal.table.Treat2_out$b_fstat_pval <- b_fstat_pval_tojoin.Treat2

for (i in 2:ncol(bal.table.Treat2_out)){
  bal.table.Treat2_out[,i] <- as.character(round(as.numeric(bal.table.Treat2_out[,i]),3))
  if (i <= 6){
    for (j in 1:nrow(bal.table.Treat2_out)){
      if (j %% 2 == 0){
        bal.table.Treat2_out[j,i] <- paste0("(",bal.table.Treat2_out[j,i],")")
      }
    }
  }
}
colnames(bal.table.Treat2_out) <- c("", "Bilingual", "English only",
                                    "Bahasa Melayu only",
                                    "F Statistic", "Prob > F")

##### RESULTS (Table 3 and Table OA3e.32) ##### 
m1 <- lm(EthnicDis ~ Treat1, data = MY_experiment)
coeftest(m1, vcov. = vcovHC(m1, type = "HC0"))

m2 <- lm(EthnicDis ~ Treat1 +
           Female + Education +
           Income + Age +
           EngFluent + MalayFluent +
           factor(Region),
         data = MY_experiment)
coeftest(m2, vcov. = vcovHC(m2, type = "HC0"))

##### DISAGGREGATING TREATMENT VARIABLE (Table OA3f.33) ##### 
###### Reference = Bilingual ###### 
m1.3_cat <- lm(EthnicDis ~ Treat2, data = MY_experiment)
coeftest(m1.3_cat, vcov. = vcovHC(m1.3_cat, type = "HC0"))

m2.3_cat <- lm(EthnicDis ~ Treat2 +
                 Female + Education +
                 Income + Age +
                 EngFluent + MalayFluent +
                 factor(Region),
               data = MY_experiment)
coeftest(m2.3_cat, vcov. = vcovHC(m2.3_cat, type = "HC0"))

###### Reference = Bahasa Melayu ###### 
MY_experiment$Treat3 <- factor(MY_experiment$Treat2, levels = c("Malay", "English", "Bilingual"))

m1.3_cat_Malay <- lm(EthnicDis ~ Treat3,data = MY_experiment)
coeftest(m1.3_cat_Malay, vcov. = vcovHC(m1.3_cat_Malay, type = "HC0"))

m2.3_cat_Malay <- lm(EthnicDis ~ Treat3 +
                       Female + Education +
                       Income + Age +
                       EngFluent + MalayFluent +
                       factor(Region),
                     data = MY_experiment)
coeftest(m2.3_cat_Malay, vcov. = vcovHC(m2.3_cat_Malay, type = "HC0"))

##### PERSPECTIVE-TAKING MECHANISM (Table OA3g.34) ##### 
m3 <- lm(PT ~ Treat1, data = MY_experiment)
coeftest(m3, vcov. = vcovHC(m3, type = "HC0"))

m4 <- lm(PT ~ Treat1 +
           Female + Education +
           Income + Age +
           EngFluent + MalayFluent +
           factor(Region),
         data = MY_experiment)
coeftest(m4, vcov. = vcovHC(m4, type = "HC0"))









